# pandas
import coptpy as copt
import pandas as pd

# copt
import coptpy as cp
from coptpy import COPT

from collections import defaultdict


# 读取csv
def read_csv(filename):
    try:
        df = pd.read_csv(filename, encoding='utf-8')
    except UnicodeDecodeError:
        # 如果 utf-8 编码失败，尝试使用其他编码
        # 或者使用 'latin1' 或 'gbk'
        df = pd.read_csv(filename, encoding='ISO-8859-1')
    return df


data_parser_folder = 'data_parser/'

# 接下来是各种图的信息，先读节点信息
df_node = read_csv(data_parser_folder + 'opt_node.csv')

# 再读边信息
df_edge = read_csv(data_parser_folder + 'opt_edge.csv')

# 最后读距离时效信息
df_order = read_csv(data_parser_folder + 'opt_order.csv')

# 开始为计算需要的各种变量赋值
# 以下变量公式中i j下标从1开始，为了代码方便，下面的i j都从0开始
# 直接已知量（可以读出来）
Dx = defaultdict(dict)  # CDC与RDC间的距离
Tx = defaultdict(dict)  # CDCi向RDCj调拨用时
Dy = defaultdict(dict)  # RDC与门店间的距离
Ty = defaultdict(dict)  # RDCi向门店j调拨用时
IT = defaultdict(dict)  # 从CDC仓库(i)到门店(j)的调拨用时是否小于10小时
M = {}  # 门店进口水果需求
N = {}  # 门店内陆水果需求
F = {}  # 仓库开设成本
G = {}  # 仓库处置成本
S = {}  # 仓库最大处置量
# 间接已知量（可以计算出来）
K1 = 0  # RDC的数量
K2 = 0  # 门店的数量

# Q2.1 的两个变量，在此修改二者的值以采样
TT = 600
RR = 0.5

# 下面的ID均指文件中的原始Id，Index均指连续编号
CDCId2Index = {}  # CDCId到Index的映射
Index2CDCId = {}  # Index到CDCId的映射

RDCId2Index = {}  # RDCId到Index的映射
Index2RDCId = {}  # Index到RDCId的映射

DCId2Index = {}  # DCId到Index的映射
Index2DCId = {}  # Index到DCId的映射

STROEId2Index = {}  # STROEId到Index的映射
Index2STROEId = {}  # Index到STROEId的映射
Id2Type = {}  # Id到Type的映射

# 未知量在后面

# 已知量赋值
# 从node获得F G S K1 K2以及映射关系
for _, row in df_node.iterrows():
    id = int(row['Id'])
    type = row['Type']
    Id2Type[id] = type
    if type == 3:
        STROEId2Index[id] = K2
        Index2STROEId[K2] = id
        K2 += 1
    elif type == 2:
        RDCId2Index[id] = K1
        Index2RDCId[K1] = id
        K1 += 1
    elif type == 0:
        # 内陆CDC
        CDCId2Index[id] = 0
        Index2CDCId[0] = id
    elif type == 1:
        # 港口CDC
        CDCId2Index[id] = 1
        Index2CDCId[1] = id

# 合并RDCId2Index和CDCId2Index得到DCId2Index
# 同时把重复的重新编号加到后面
DCId2Index = {**RDCId2Index,
              Index2CDCId[0]: K1, Index2CDCId[1]: K1+1}
# 合并Index2CDCId和Index2RDCId得到Index2DCId
Index2DCId = {**Index2RDCId,
              K1: Index2CDCId[0], K1+1: Index2CDCId[1]}

#重新回炉一遍node把F G S对应好
for _, row in df_node.iterrows():
    id = int(row['Id'])
    type = row['Type']
    # 填充F G S
    if type == 3:
        continue;

    F[DCId2Index[id]] = row['Opening_fee']
    G[DCId2Index[id]] = row['Processing_fee']
    S[DCId2Index[id]] = row['Capacity']

#输出DCId2Index全部内容
print(DCId2Index)
print(Id2Type)

sum_M = 0
sum_N = 0

# 从order获得M和N
for _, row in df_order.iterrows():
    id = row['Id']
    N[STROEId2Index[id]] = row['dm_qty']
    M[STROEId2Index[id]] = row['im_qty']

    sum_M += row['im_qty']
    sum_N += row['dm_qty']
    
# 从edge获得Dx Tx Dy Ty
for _, row in df_edge.iterrows():
    fromId = row['From']
    toId = row['To']

    # 根据类型决定填充谁，有CDC向RDC和CDC向门店和RDC向门店三种情况
    if Id2Type[fromId] <= 1 and Id2Type[toId] == 2:
        # CDC向RDC
        Dx[CDCId2Index[fromId]][RDCId2Index[toId]] = row['Distance']
        Tx[CDCId2Index[fromId]][RDCId2Index[toId]] = row['Duration']
    elif Id2Type[fromId] <= 1 and Id2Type[toId] == 3:
        # CDC向门店
        Dy[DCId2Index[fromId]][STROEId2Index[toId]] = row['Distance']
        Ty[DCId2Index[fromId]][STROEId2Index[toId]] = row['Duration']
    elif Id2Type[fromId] == 2 and Id2Type[toId] == 3:
        # RDC向门店
        Dy[DCId2Index[fromId]][STROEId2Index[toId]] = row['Distance']
        Ty[DCId2Index[fromId]][STROEId2Index[toId]] = row['Duration']

#赋值IT
for DCIndex in Index2DCId:
    for STOREIndex in Index2STROEId:
        if(Ty[DCIndex][STOREIndex] <= TT):
            IT[DCIndex][STOREIndex] = 1
        else:
            IT[DCIndex][STOREIndex] = 0


# 创建未知量，主要包括：
# Xij ： CDCi向RDCj调拨量
# Yij： RDC+CDC i向门店j调拨量
# Hi : RDC+CDC i是否智能采购
# Li ： RDC i是否存在
# 创建模型# 初始化COPT环境和模型
# Create COPT environment
env = cp.Envr()
model = env.createModel("lp_ex1")

model.setLogFile("log.txt")
model.setParam(copt.COPT.Param.GPUMode, 0)
model.setParam(copt.COPT.Param.LpMethod, 6)
model.setParam(copt.COPT.Param.RelGap, 0.05)
model.setParam(copt.COPT.Param.AbsGap, 1.0)
model.setParam(copt.COPT.Param.TimeLimit, 5*60.0)  # 限制求解时间为300秒，求一个粗略解即可


# 定义变量
x = model.addVars(2, K1, vtype=COPT.CONTINUOUS, nameprefix='x')
y = model.addVars(K1+2, K2, vtype=COPT.CONTINUOUS, nameprefix='y')
l = model.addVars(K1, vtype=COPT.BINARY, nameprefix='l')
H = model.addVars(K1+2, vtype=COPT.BINARY, nameprefix='H')
h = model.addVars(K1+2, vtype=COPT.BINARY, nameprefix='h')
q = model.addVars(K1+2, K2, vtype=COPT.BINARY, nameprefix='q')
r = model.addVars(K1+2, K2, vtype=COPT.BINARY, nameprefix='r')
Q = model.addVars(K1+2, K2, vtype=COPT.BINARY, nameprefix='Q')
SS = model.addVars(K1, vtype=COPT.CONTINUOUS, nameprefix='SS')
p = model.addVars(K1+2, vtype=COPT.CONTINUOUS, nameprefix='p')
Q_y = model.addVars(K1+2, K2, vtype=COPT.CONTINUOUS, nameprefix='Q_y')
T_y = model.addVars(K1+2, K2, vtype=COPT.CONTINUOUS, nameprefix='T_y')
Y = model.addVars(K1+2, K2, vtype=COPT.CONTINUOUS, nameprefix='Y')

# 定义目标函数
term1=cp.quicksum(0.6 * x[i, j] * Dx[i][j] for i in range(2) for j in range(K1))
term2=cp.quicksum(1.25 * y[j, k] * Dy[j][k] for j in range(K1+2) for k in range(K2))
term3=cp.quicksum(F[i] * l[i] for i in range(K1))
term4=cp.quicksum(0.5*G[i]*H[i]+10000*h[i] for i in range(K1+2))
term5=cp.quicksum((p[i] - H[i])*G[i] for i in range(K1+2))

objective = term1 + term2 + term3 + term4 + term5

model.setObjective(objective, COPT.MINIMIZE)

# 添加约束条件# 添加约束条件
BIG_M = 1e8  # 极大常数

# Q_{ij} \le q_{ij}
model.addConstrs((Q[i, k] <= q[i, k] for i in range(K1+2)
                 for k in range(K2)))

# Q_{ij} \le r_{yij}
model.addConstrs((Q[i, j] <= r[i, j] for i in range(K1+2)
                 for j in range(K2)))

# Q_{ij} \ge q_{ij} + r_{yij} - 1
model.addConstrs((Q[i, j] >= q[i, j] + r[i, j] - 1 for i in range(K1+2)
                 for j in range(K2)))

model.addConstrs((cp.quicksum(Q[i, k] for i in range(K1+2))<= 1 for k in range(K2)))

model.addConstrs((cp.quicksum(q[i, k] - Q[i, k] for i in range(K1+2)) <= 1 for k in range(K2)))

model.addConstrs((cp.quicksum(q[i, j] for i in range(K1,K1+2)) == l[j]
                 for j in range(K1)))

# x_{ij}, y_{ij} \ge 0
model.addConstrs((x[i, j] >= 0 for i in range(2)
                 for j in range(K1)))
model.addConstrs((y[j, k] >= 0 for j in range(K1)
                 for k in range(K2)))
biggest_MN = sum_M+sum_N


# x_{ij}, y_{ij} \ge 0
model.addConstrs((x[i, j] <= biggest_MN for i in range(2)
                 for j in range(K1)))
model.addConstrs((y[j, k] <= biggest_MN for j in range(K1)
                 for k in range(K2)))

# \sum_{i = 0}^{1}x_{ij} = \sum_{j = 0}^{K_1-1}y_{ij}
model.addConstrs((cp.quicksum(x[i, j] for i in range(2)) == cp.quicksum(
    y[j, k] for k in range(K2)))for j in range(K1))

# p_i = \sum_{j = 0}^{K_2-1}y_{ij}
model.addConstrs((p[j] == cp.quicksum(y[j, k] for k in range(K2))
                 for j in range(K1)))

# SS_i = s_il_i

# S_i \ge s_i - (1 - l_i)M
model.addConstrs((SS[j] >= S[j] - (1 - l[j]) *
                 BIG_M for j in range(K1)))

# S_i \le s_i
model.addConstrs((SS[j] <= S[j] for j in range(K1)))

# S_i \le l_iM
model.addConstrs((SS[j] <= l[j] * BIG_M for j in range(K1)))

# S_i \ge -l_iM
model.addConstrs((SS[j] >= -l[j] * BIG_M for j in range(K1)))

# p_i \le S_i
model.addConstrs((p[j] <= SS[j] for j in range(K1)))

# y_{ij} \le Mq_{ij}
model.addConstrs((y[j, k] <= BIG_M * q[j, k] for j in range(K1)
                 for k in range(K2)))

# Q_{yij} = Q_{ij}y_{ij}

# Q_{yij} \ge y_{ij} - (1 - Q_{ij})M
model.addConstrs((Q_y[j, k] >= y[j, k] - (1 - Q[j, k]) *
                 BIG_M for j in range(K1+2) for k in range(K2)))

# Q_{yij} \le y_{ij}
model.addConstrs((Q_y[j, k] <= y[j, k] for j in range(K1+2)
                 for k in range(K2)))

# Q_{yij} \le Q_{ij}M
model.addConstrs((Q_y[j, k] <= Q[j, k] * BIG_M for j in range(K1+2)
                 for k in range(K2)))

# Q_{yij} \ge -Q_{ij}M
model.addConstrs((Q_y[j, k] >= -Q[j, k] * BIG_M for j in range(K1+2)
                 for k in range(K2)))

model.addConstrs((cp.quicksum(Q_y[i, k] for i in range(K1+2)) >= M[k] for k in range(K2)))

# H[j]系列的限制
model.addConstrs(H[j]>=p[j]-(1-h[j])*BIG_M for j in range(K1 + 2))

model.addConstrs(H[j]<=p[j] for j in range(K1 + 2))

model.addConstrs(H[j]<=h[j]*BIG_M for j in range(K1 + 2))

model.addConstrs(H[j]>=-h[j]*BIG_M for j in range(K1 + 2))

# Y_{ij} = q_{ij}y_{ij}

# Y_{ij} \ge y_{ij} - (1 - q_{ij})M
model.addConstrs((Y[j, k] >= y[j, k] - (1 - q[j, k]) *
                 BIG_M for j in range(K1+2) for k in range(K2)))

# Y_{ij} \le y_{ij}
model.addConstrs((Y[j, k] <= y[j, k] for j in range(K1+2)
                 for k in range(K2)))

# Y_{ij} \le q_{ij}M
model.addConstrs((Y[j, k] <= q[j, k] * BIG_M for j in range(K1+2)
                 for k in range(K2)))

# Y_{ij} \ge -q_{ij}M
model.addConstrs((Y[j, k] >= -q[j, k] * BIG_M for j in range(K1+2)
                 for k in range(K2)))

# \sum_{i = 0}^{K_1 + 1}(Y_{ij} - Q_{yij}) \ge n_j
model.addConstrs((cp.quicksum(Y[i, k] - Q_y[i, k] for i in range(K1+2)) >= N[k]
                 for k in range(K2)))

# T_{ij} = 1 if t_{yij} \le 10 else 0
# T被IT替代

# T_{yij} = T_{ij}y_{ij}

# T_{yij} \ge y_{ij} - (1 - T_{ij})M
model.addConstrs((T_y[j,k] >= y[j, k] - (1 - IT[j][k]) *
                 BIG_M for j in range(K1+2) for k in range(K2)))

# T_{yij} \le y_{ij}
model.addConstrs((T_y[j,k] <= y[j, k] for j in range(K1+2)
                 for k in range(K2)))

#下面两个是常量约束常量，可以去掉                 
# T_{yij} \le T_{ij}M
model.addConstrs((T_y[j,k] <= IT[j][k] * BIG_M for j in range(K1+2)
                for k in range(K2)))

# T_{yij} \ge -T_{ij}M
model.addConstrs((T_y[j,k] >= -IT[j][k] * BIG_M for j in range(K1+2)
                 for k in range(K2)))

# \frac{\sum_{i = 0}^{K_2}\sum_{j = 0}^{K_2}T_{yij}}{\sum_{j = 0}^{K_2}(m_j + n_j)} \ge R
f_up = cp.quicksum(T_y[i,k] for i in range(K1+2) for k in range(K2))
f_down = cp.quicksum(M[k] + N[k] for k in range(K2))
model.addConstr(f_up >= RR*f_down)

# 求解模型
model.solve()



# 输出结果
if model.status == COPT.OPTIMAL:
    print('Optimal solution found:')
    for i in range(2):
        for j in range(K1):
            print(f'x[{i},{j}] = {x[i, j].x}')
    for i in range(K1+2):
        for j in range(K2):
            print(f'y[{i},{j}] = {y[i, j].x}')
else:
    print('No optimal solution found.')

model.write("lp_ex1.mps")
model.write("lp_ex1.bas")
model.write("lp_ex1.sol")
model.write("lp_ex1.par")
